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In this article, a fast algorithm is developed to compute two-dimensional convolutions 
of an array of d } • d 2 complex number points, where d 2 = 2 m and d ; = 2 m ~ r+ 1 for some 
1 < r < m. This new algorithm requires fewer multiplications and about the same number 
of additions as the conventional FFT method for computing the two-dimensional 
convolution. It also has the advantage that the operation of transposing the matrix of 
data can be avoided. 


I. Introduction 

Two-dimensional convolutions of two sequences of complex number points can be applied to many areas, in particular to the 
synthetic aperture radar (SAR) (Refs. 1,2). In SAR, a two-dimensional cross correlation of the raw echo data of complex 
numbers with the response function of a point target is required to produce images. When the two-dimensional filter does not 
change rapidly with the range, one can divide the entire range of echo data into several subintervals. Within each subinterval, one 
can use a constant filter function. This is accomplished usually by using the conventional fast Fourier transform (FFT). However, 
the FFT algorithm generally requires a large number of floating-point complex additions and multiplications. Also, the transpose 
of a matrix is usually required in the computation of such a two-dimensional convolution. 

Recently, Rader (Ref. 3) proposed that a number-theoretic transform (NTT) could be used to accomplish two-dimensional 
filtering. It was shown (Ref. 4) that an improvement, both in accuracy and speed, of two-dimensional convolutions could be 
achieved by transforms over a finite field GF{q), where q is a prime of the form 45 X 2 29 + 1. However, to compute a 
two-dimensional convolution of two long sequences of an integer number of points, such a transform over a finite field did not 
allow for a wide variety of dynamic ranges. 
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More recently, Nussbaumer and Quandalle (Ref. 5) showed that a type of polynomial transform over the complex numbers 
could be used to efficiently compute two-dimensional convolutions. A principal advantage of this method over the above 
mentioned techniques is that the need for computing the transpose of matrix can be avoided. Furthermore, this new method 
offers a wider variety of dynamic ranges. It was shown recently by Arambepola and Rayner (Ref. 6) that the ideas of Nussbaumer 
and Quandalle can be generalized to a radix-2 polynomial transform analogous to the conventional radix-2 FET. But they do not 
make use of the Chinese remainder theorem as it is shown in this article to further reduce the complexity of the algorithm. 

In this article, it is shown that a combination of a fast polynomial transform (FPT) and the Chinese remainder theorem (CRT) 
can be used to very efficiently compute a two-dimensional convolution of a d l • d 2 complex number array, where d 2 ,= 2 m and 
c/j = 2 m ~ r+l for 1 <r<m. Such a new algorithm requires considerably fewer multiplications and about the same number of 
additions as the conventional algorithm for the two-dimensional case. Therefore, it has the potential for important application in 
SAR. 


II. The Computation of Two-Dimensional Convolutions 

The following algorithm for a two-dimensional digital convolution is based on an important identity. Let d. be a power of 2 for 
i=l, 2 and let C be the field of complex numbers. Also let a tl , t2 and b fjj(2 be two d^d 2 arrays, where O^t^d.- 1 for 
1 = 1,2. Then the two-dimensional cyclic convolution of a tl , 2 and b where a tlit2 , b ty j 2 e c, is defined by 


d r-i d 2~ l 


C n v n 2 = E E a t v t 2 b (ni -t lU n 2 -t 2 )’ 0 < n. < d. for i = 1 , 2 


*t“0 f 2 = 0 


where (n. - t.) denotes the residue of n. - t. modulo d. for /' = 1 , 2. 
Define 


a 2 -i 


d 2~ 1 


A tSV = E V,, BAZ) = 2 K.Z-- 


t 2 = o *»■ 2 


1 r 2 =o 1- 2 


d 2 -t 


C */ Z) = E C «,-«, Z2 ’ for 0<n v t x <d r 1 

* — n 12 


« 2 = 0 


0<n 2 , t 2 <d 2 - 1 


(1) 


( 2 ) 


Nussbaumer and Quandalle (Ref. 5) expressed the two-dimensional convolution in Eq. (1) as a one-dimensional convolution of 
polynomials, i.e., 



d i-i 

= E A t W B (n _ r } (Z) mod (z* 2 - l) for 0<« 1 <d l - 1 


where C ni (Z),A tl (Z),B tl (Z) are defined in Eq. (2) and (« 1 - denotes the residue of - t x modulo 


( 3 ) 
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If d 2 = 2 m and d = 2 m r+ 1 for 1 < r < m, then one can factor Z dl - 1 into pairwise relatively prime factors as follows, 

1 - (zV + 1) l)... (/•"" + 1 )(/»"'+ l)^'- 1 j 

Thus, Eq. (3) is equivalent to 


C »,® ■ 5] ^ 1 ®fi<„ i ., i) (Z)mod(z' ij,I -u)(z v2 + l) l) (4) 


Since (Z^ 2 ^ 2 + 1), ■ • • ,(Z d 2/ 2 - 1) are pairwise relatively prime, by the Chinese remainder theorem (CRT) for polynomials 
(Ref. 7), the polynomial congruences 


C n ( Z) = C n (Z) mod ^Z* 2 ' 2 + lj fori =1,2, ••■,/■ 
and 

C «j(Z) = C (Z) mod (z^ 2 - l) 


(5a) 


(5b) 


have a unique solution C n ^ (Z) given by 


Z^ 2 


+ C (Z) 


z d ^ 2 ' + 1 


(v) 


z** 2 - 1 
Z W - 1 


mod Z 2 - 1 


( 6 ) 


Note that the arithmetic needed to compute Eq. (6) requires only cyclic shifts and additions. The number of real additions needed 
to compute Eq. (6) is 2 rd^d^. 


The derivation of C‘ t (Z) in Eq. (5a) proceeds in the following manner: 


C' (2) = C (Z) mod 
l l 



/ d -> i 2 ‘ \ 

3 T, A ,w\»,-d z >™ i [ z +1 ) 

fj-0 


i_> / d J2 1 \ 

s E a ; mt,, -< 1 >< z)mod z +1 < 7 > 

r t =o iii \ ' 
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where A ' (Z) is defined by 
1 


(d 2 /2')-l 


A h (Z) ,? 0 ( V2 " VW 2 * + % *2 +2< V 2 ' 


% f 2 +3< V 2,+ ' 




(«yn-i 


= X! a' Z 2 for l<i<r 

l' r 2 


V° 


and B\^(Z) is obtained from the expression /I '^(Z) on replacing each a. . by b . that is, 

(d 2 / 2')-i 

5' (Z) = Y* hi , z' 2 for l</<r 

f l *1**2 

* 2 = ° 1 2 

Note that A\^{Z) and B l f (Z) are, respectively, A t ^(Z) and5^(Z) reduced modulo(Z d2 ^ 2 + 1). 

It was demonstrated in (Ref. 6) that a fast polynomial transform can be used to compute Eq. (7). Let us show this more 
carefully. Since d 2 = 2 m and d 1 = 2 m_r+1 for 1 < r <m, then dj 2 i_1 = 2 m_,+1 = 2 r ~‘ • 2 m_r+1 = 2 r ~ i d l for i = 1, 2, ■ • • , r. 
Note the identity 


z' 2 * s = (-1) 2 — 1 mod (z?* 12 +lj 

Also let / = 2 r ~‘ and define the d -point / th power polynomial transform ofA‘ (Z) and B! (Z) by 

i r i r 2 

d j-l t (d 2 /2 f )-l , . 

A^(Z) s £ A' f (Z)(Z) fl * = X] aj^rnod fZ* 2 * +1 forO<* <<*,- 1 

t =0 1 fi=0 ’ ' 


and 


V 1 (d 2 ft‘)-i i 

B[ (Z) s £ B\ (Z)(z/ lfc = 2 ^' M ZSmod r 2 +1 ) f°rO 1, 

r,=o 1 e=o ’ ' ’ 


( 8 ) 


(9a) 


(9b) 


respectively, where Z = Z 1 is the ; th power of Z. Note by Eq. (8) that 


~d, id , 2 r 'd, d /2* 1 / d_/ 2 ' \ 

Z 1 = Z 1 =Z 1 =Z 2 = 1 mod (z 2 +11, 


for / = 1, 2, • • • ,r 
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The product of the transforms A‘ k (Z) and B‘ k (Z) is given by 


(<y2')-t 


(Z) = A‘ k (Z) • B‘ k (Z) = £ ? k ' S Z 2 mod (z^ + l) 

e=o ' ' 


One needs now to compute the inverse transform of C'(Z). That is, to compute 


forO <k<d 1 - 1 


1 <i<r 


V 1 




1 k= 0 


( 10 ) 


d r 1 l d \~ x 


< d i - 1 


= ^i 1 DIE A i& zsk I • ( E Kw * ) 2~*" 1 

fc =0 \ i =0 / \ r =0 


= E E ^*,(2)^(2) J- E Z^ r - n ^ mod Iz^ 2 ' 2 +l) 

i=0 r= 0 1 fc=0 ' ' 


Now let t =s + r- n 1 and note that 


( 11 ) 


~fd. 

1 Z 1 - 1 


o = _L v jtk _ i z 1 - i 

1 (S') - 1 _ 1 (z d ^ ) - 1 _ / d,/2'' \ 

?-i = ^7 z* - \ — = 0mod r +1 ) for ^ 0modd i 

It is seen that 5 = 1 for t = 0 mod d y and S = 0 for t ^0 mod d l . Hence Eq. (1 1) yields the desired result, namely, 


dj-l . 

<2,® a £ 4 1 < z)B '<»,-v <z>n ’° d ( Z<V * >) 


forO<rZj <d 1 - 1, 1 </</• 


( 12 ) 


where (n 1 - r ( ) denotes the residue of(/jj - fj) modulodj. 

Since d 2 = 2 m and d y = 2 m ~ r+1 for 1 < r^m the analogue of the conventional FFT algorithm can be utilized to realize the 
needed polynomial transforms Eqs. (9a) and (9b) and its inverse defined by Eq. (11). That is, in Eq. (9a) let the input data to the 
usual FFT be replaced by a sequence of polynomials over C, and let w, the usual d^-th complex root of unity, be replaced by Z, 
satisfying Z^ 1 = 1, and finally let the arithmetic operations on the complex number field be replaced by arithmetic operations in 
the field of polynomials. Then the polynomial transforms in Eqs. (9a), (9b), and (1 1) can be computed by the conventional FFT 
algorithm. The polynomial transform, obtained in this manner, is called the fast polynomial transform (FPT). The number of 
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operations needed to perform the polynomial transform defined in Eq. (9a) requires only (d 2 /2) ■ log2, cyclic shifts of 
polynomials and 2(d 2 l2)d 1 -\ogd l real additions. 

To compute C k (Z), given in Eq. (10) for 0 - 1, one needs to compute a polynomial product modZ^ 2 ^ 2 + 1. To do 

this directly takes a lengthy calculation. Using an idea, based on the work of Arambepola and Rayner (Ref. 6), this equation can 
be transformed into 


C k ( cu ) = A‘ k (cu) • B‘ k ( cu ) mod 



(13) 


d 12* 

by changing variables from Z to u with the mapping Z = cu, where c is a dj2-th primitive root of - 1, i.e. c 2 = - 1. One can 

verify (Ref. 8) that the above polynomial product mod u dl ^ 2 ‘ - 1 can be obtained as a cyclic convolution of the twoc/ 2 /2'-point 
coefficients of the polynomials A k (cu) and B‘ k (cu). This is vastly simpler than finding the product modulo Z d 2 / 2 ' + 1. 


In Eq. (13), one observes that 


(rf 2 /2')-l 


( V 2 ')-1 


A‘ k (cu) = Y (cuf = Y, K,i u * 


8 = 0 


(d 2 /2‘)-l 


8=0 


(d 2 /2')-l 

B^icu) = Y **.«<«)“ = E %,s 


8=0 


8=0 


(14) 


and 


(d 2 /2‘)-l 


(d 2 l 2‘)-l 

C‘(cu) = Y M 2 = E c‘ k 8 « e 


8=0 


8=0 


where 2^ = cV fcE , t‘ ki = and c‘ fc g = c s c' fc s for 0 < £ < (dj 2') - 1. If one defines 


(d 2 /2')-l 


(<V 2 ')-i 

A‘ k ( u ) = Y K («) = E K.t 


8=0 


8=0 


(rf_/2')-l 


Cl(u) = Y 


8=0 


(15) 


Then Eq. (13) becomes 
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C‘ k (u) = A l k {u) B‘ k (u) mod i^u 2,2 - lj 


(16) 



The polynomial multiplication in Eq. (16) can be computed by the cyclic convolution of the two £? 2 /2'-point sequences, 
‘ >“k,d 2 i 2 '-i) and (&'*.„, , • • • ,b i ktd2 / 2 ‘- 1 )- This cyclic convolution can be written as 


(<y 2 ')-i 

E for ° 


0<Z<d 2 /2 i - 1 


where (£ - n) denotes the residue of (£ - n) modulo d 2 l2‘ and a‘ n k and b‘ n k are the coefficients of the A‘ k (Z) and B' k (Z), 
respectively. The convolution c' k i in Eq, (17) can be computed by using the conventional FFT (Ref. 9). Thus, the desired 

_ • r? / -t r\ • * 'Vi' . a ^ f. v j i - n - /a/ i ____ _ 1 _a * 1 f ^ 'V — Q-O-i »t*i. _ -C „__1 


The flow chart for computing C*(Z) is given in Fig. 1. 
Consider now the computation of C*^(Z) in Eq. (5b). That is, 


d i 1 / d /2 r \ 

C® ■ E < (Z) > (z)modz! -> 

1 f =0 1 11 ' ' 


(d 2 /2 r )-l 

^ f i (Z)= f ? 0 K> f 2 +af r f 2 +<, 2/ 2,+flf iW 


+ a + . . . + a 

t v t 2 +3d 2 i2 r t v t 2 H2 r -i)d 2 /2 r 


(d /2 r )-l 

- e *;,, z2 

' *,,< 2 

f 2 =0 

and is obtained from the expression A^(Z) on replacing each a t j by by, that is, 


(d 2 /2 r )-l 


B\ (Z) = £ K t z '° 

i rri 2 


Evidently, C* (Z) in Eq. (18) is in the form of 


(d 2 /2 r )-l 


c' n (z) = X) < „ z ' 2 

"l __ "l ,n 2 

n 2 = 0 

* . rf,/2' 

Note here that f (Z)and2? f (Z) are, respectively, ,4^ (Z) and (Z) reduced modulo Z + 1. 
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In order to use the FPT technique to obtain Eq. (18), again use the idea of Arambepola and Rayner (Ref. 6). To compute Eq. 
(18), letZ = cu, where c is a d 2 l2 r -th root of - 1. By an argument similar to that used to obtain Eq. (12), (18) becomes 

dji r \ 

; 2 + l) for 0 - 1 (19) 

In Eq. (19), one observes that 


V 1 / 

</») = Z mod ( 


and 


a; (cu) = 

i 

(d 2 /2 r )-l 

Z \, f 2 ^ 2 = 

t 2 = 0 

(d 2 /2 r )-l 

V 1 ' *2 

a t t u 

t 2 =o 


(<y2 r )-i 

(<y 2 r )-i 

B" (cu) = 

r i 

Z b 't v t 2 W 2 = 

t 2 = 0 

Z K t 
1,2 


(d 2 /2 r )-l 

< d 2/ 2 '')- 1 

c; (cu) = 

"i 

E v 2 m’ = 
" 2 =0 

Y c ' u 

^ n v n 

n 2 = 0 

t 2 * 

. =c *a t 
2 *!• 

( , and b\ , -c 2 b* t t 
* 2 f l> f 2 f l* f 2 

for0< nj,t, <dj - 


( 20 ) 


In Eq. (20), define A' ti (u) -Af^cu),^ (u) = B* t ^(cu), and C^(u) = C* ^(cw). Then Eq. (19) becomes 


C 'n = Z A 't - t )(«)mod in 2 ' 1 + lj for 0 <«j - 1 

1 fj=o 111 ' ' 

Using the same procedure used in the computation of Eq. (7), one obtains C' n ^ (u). Hence, e* in Eq. (18) is obtained from 
the substitutions, c n2 c' n n . The number of real multiplications and real additions needed to compute Eq. (19) are d x d 2 j2 r ~ 2 + 
d 1 d 2 l2 r ~ 2 log (d 2 /2 r ) + 16 d t (d 2 l 2 r - 1) and 2 d 2 d 2 + d x (d 2 f 2 r ~ 2 ) log d x + 3 d^dj 2 r ~ 1 ) log (d 2 l 2 r ), respectively. Hence the 
total number of real multiplications and real additions needed to compute Eq. (4) are 8 d { [(d 2 /2 r+1 )(2 r (log d 2 +3) +4) - 
(d 2 + r + 2)] and 2 • d l [d 2 (2r-A) + d 2 ( 2 • log <7, - 1 + 3 • log d 2 ) + 6 • d 2 /2 r ] , respectively. The flowchart of this new algorithm 
is shown in Fig. 2. 

In the introduction it was stated that transposition of the data matrix, usually required in two-dimensional convolution, can be 
avoided in this new algorithm. Assume the data matrix for a typical array such as a, ^ , 0 < t i < d t - 1 for i = 1 ,2, is arranged 

with t x indexing the row, t 2 indexing the column. Then the polynomial A t (Z) of Eq. (2) is the t x th row of this matrix. The FPT 
of A t ^(Z) can be implemented by a decimation-in-time algorithm analogous to the well-known conventional FFT one. An 
example is given in the Appendix. This requires that only two polynomials (or rows) need be available for processing at anytime 
and the resulting two polynomials can replace the input ones. Consequently, no additional storage is required for the FPT, and the 
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data array need only be accessed by rows with replacement after processing, thus obviating the accessing of individual columns 
usually done in two dimensional convolution. 


In the Appendix, it is shown by an example how fast polynomial transforms can be combined with FFTs to yield a new fast 
algorithm for computing a two-dimensional convolution. The number of operations needed for this new algorithm to perform the 
two-dimensional convolutions of ad 1 -d 2 array, where d 2 = 2 m and d t =2 m ~ r+l for 1 <r<m, is given in Table 1. In this table, 
the FPT-FFT algorithm and conventional FFT algorithm for computing the two-dimensional convolutions are compared by giving 
the number of operations to perform these algorithms. 


Table 1. Complexity of new algorithm for two-dimensional convolutions 


Convolution size 
d = dy X c? 2 = 2 m ~ r+l X 2 m 
for 1 < r < m 


FPT-FFT-algorithm 


Radix-2 FFT algorithm 

No. of factors of 
z d 2-l=r+l 

No. of real 
multiplication 

No. of real 
addition 

No. of real multiplication 
4djd 2 (log d 2 + 
logdj + 1) 

No. of real addition 
6dyd 2 (log d 2 
+ logcf j) + 2 dyd 2 

2 7 X 2 9 

4 

2,747,392 

5,537,792 

4,456,448 

6,422,528 

2 7 X 2 10 

5 

5,892,096 

12,288,000 

9,437,184 

13,631,488 

211x2" 

2 

234,831,872 

457,179,136 

385,875,968 

562,036,736 

2 10 X 2 U 

3 

109,019,136 

222,298,112 

184,549,376 

268,435,456 

2 8 X 2 12 

6 

55,035,904 

118,882,304 

88,080,384 

127,926,272 

2 8 X 2 13 

7 

117,948,416 

258,342,912 

184,549,376 

268,435,456 

2 7 X 2 13 

8 

58,842,112 

129,073,152 

88,080,384' 

127,926,272 


95 







Fig. 1 . Computation of the cyclic convolution of two Z-polynomials 
mod Z? 2 '^ + 1, i.e., (Z) given in Fig. 2 by 

fast polynomial transforms 
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Appendix 

Example of a Two-Dimensional Convolution of 
Dimension Computed by the New Method 


Example: Compute the two-dimensional cyclic convolution of two sequences ^ and b t v t 2 ' i.e., 


4-1 


4-1 


c = y y a. . b, ... . , forO<«,,n, <3 

n l' n 2 “ 1 f l’ f 2 ( "l f 2^ 1 2 


(A-l) 


* 1=0 f 2 =° 


where (x) denotes the residue of x modulo 4 and where a Q0 =y 00 = 1, a 01 = b 0>1 = l,a 0 2 = b 0 2 = 0, a Q3 = b 0 3 = 0, 

*1,0 =*1,0 =1 > fl l,l =*1.1 = 1 ’ a i,2 =b U2 = ^l.S^l.S =°>*2,0 = *2,0 =0 >* 2>1 = *2,1 = °> *2,2 = *2,2 = ^’*2,3 = 
*2,3 = 0 >*3,O = *3,0 =0 > *3,1 = *3,1 =0 >*3,2 = *3,2 = °> *3,3 = *3,3 =0 - 


Define 


f,=o 1 f =0 1 2 


‘2 

4-1 


C„,(Z) = for0<»,^,<3 

1 « -n 1 * 


n 2 =0 


That is,^ 0 (Z) = J B 0 (Z) = 1 + Z,Aj(Z) =I?j(Z) = 1 + Z + Z 2 ,A 2 (Z) =B 2 (Z) = 0,A 3 (Z) =B 3 (Z) = 0. From Eq. (4), (A-l) becomes 


c n V) = 


4-1 

S A, (Z) B (n } (Z) mod (Z 2 - 1)(Z 2 + 1) 

* -n t 11 


From Eqs. (5a) and (5b) one obtains 


4-1 

C‘ (Z) = Y, A] (Z) ^ _ f } (Z) mod (Z 2 + 1) 

1 + -A 1 


where A l Q (Z) = Z^(Z) = 1 +Z,ylj(Z) = i?j(Z) = Z,A 2 (Z) = 5»(Z) = 0, A \{Z) = Z?‘(Z) = 0. 
and 


(A-2) 


(A-3a) 


4-1 

C„‘(Z) = < Z > 5 (‘„ -r ) (2) mod ( z2 - !) 

r,=o 1 1 1 


(A-3b) 


where A* 0 (Z) = B* 0 (Z) = 1 + Z, A\{Z) = B\(Z) = 2 + Z,^ *(Z) = 5* (Z) = 0, A;(Z) = 5^(Z) = 0. 
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To compute Eq. (A-3a), one takes the Z polynomial transform of 4? (Z) given by 

r i 

A k (Z) .= £ mod ( Z2 + 1) 


The flowgraph of the decimation-in-time decomposition of a four-point polynomial transform computation is given in Fig. 1-A as 
A 0 (Z) = 1 + 2 Z, A j (Z) = Z,A 2 (Z) = 1, and A 3 (Z) = 2 + Z. Similarly, the Z polynomial transform of B k (Z) are B 0 (Z) = 1 + 2Z, 
B , (Z) = Z, B 2 (Z) = 1 , and B 3 (Z) = 2 + Z. 

Let 


C k (Z) = A k (Z) • B k (Z) mod (Z 2 + 1) for k = 0, 1 , 2, 3 (A-4) 

If one uses the mapping Z = iii, where i 2 = - 1 then Eq. (A-4) becomes 

C k (iu ) = A k (iu) ■ B k (iu ) mod (u 2 - 1) for 0 < k < 3 (A-5) 

For t = 0, Eq. (A-5) becomes 

C 0 (iu) = A Q (iu) ■ B Q (iu) = (1 + 2/«)(l + 2 hi) mod (« 2 - 1) (A-6) 


The FFT can be used to compute Eq. (A-6). That is, let a 0 = l,a l =2i,b 0 = l,b 1 = 2i. The transform of a n is 


A 


k 


2-1 nk 2-1 

H a n W = Y. a n^ nk=l+ ^ k f0r * = °. 1 

k= 0 n=0 


(A-7) 


The radix-2 FFT algorithm is used to compute Eq. (A-7). The results are A 0 = 1+2 i and A x = 1-2 i. Similarly, one obtains B 0 = 
1 + 2/ and5 t = 1-2 i. Let C 0 =A 0 B 0 = - 3 + 4 i, and =A 1 •B 1 =-3 - 4i. The inverse Fourier transform of C k is 




k= 0 


= 2~ 1 




k=0 


= 2 -1 [(-3 + 4z) + (-3 - 40 (- 1)“" ] for n = 0, 1 


(A-8) 


Again, the radix-2 FFT algorithm is used to compute Eq. (A-8). That is, c 0 = -3 and c v = 4 i. Thus, C 0 (iu ) = - 3 + Aiu mod (u 2 - 
1). Hence, one obtains C 0 (Z) = -3 + 4/ -, rZ = -3 + 4Z mod (Z 2 + 1). 

Similarly one obtains C,(Z) = - 1, C 2 (Z) = 1, and C 3 (Z) = 3 + 4Z. The inverse polynomial transform of C k (Z) is given by 


C' n (Z) = 4 -1 

"t 


E ? fc (Z) Z _ " 


(A-9) 


The FPT is used to compute Eq. (A-9). The results in Eq. (A-3) are C*(Z) = 2Z, Cj(Z) = - 2 + 2Z, C 2 (Z) = - 1 , and C*(Z) = 0. 
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To compute Eq. (A-3b), one needs to use the mapping Z = iu. Then Eq. (A-3b) becomes 


4-1 



(n 1 -t 1 )0“) 


= c* +ic* u mod(« 2 + l) for0<n.<3 
"i>o ”i-i 1 


(A- 10) 


where A *(/'«) = 2?*(/u) = 1 + iu, A | (iu) = B*(iu ) = 2 + iu, A* 2 (iu ) = /?*(/«) = 0,yf*(;«) = = 0. 


Let (u) =B'(u) = 1 + iu, A ' (w) = B \ (») = 2 + iu, A'Ju) = B' 2 (u) = 0,A 3 (u) = B'Ju) = 0 and Cf n (u) = c n +c'„ u for 0 < 

o ' u * j > _ _• * it /a inn 1 


n, < 3, where c„ 

1 7 n \ 


= c„ 


' 1,0 


and c 


’ 1,1 


= ic 


" 1,0 


. Hence, Eq. (A-10) becomes 


4-1 

C n = X) -r )( u ) m °d (m 2 + 1) (A-ll) 

1 ij=o 1 1 1 

Equation (A-l 1) can be computed by the FPT. This yields Cf Q (u) = 2 + 2 iu, C^iu) = 6 + 6 iu, C? 2 (u) = 5 + 4 iu, and C^n) = 0. 
Hence 


C*(Z) = d 0 (C'Z) = 2 + 2 Z, C*(Z) = C[(T l Z) = 6 + 6 Z ; 

c*(z) = q'o-'z) = 5 + 2Z, c*(z) = c^r'z) = o (A- 12 ) 

From Eq. (6), one obtains 

C n &) = \ \d ^ (2) (Z 2 + 1) + C\(Z) (1 - Z 2 )] mod (Z 4 - 1) 

for 0 <^<3 


where C t ] (Z) and C* (Z) are defined in Eqs. (A-3a) and (A-3b), respectively. Thus, the desired results are C Q (Z) = 1 + 2Z + Z 2 , 
C t (Z) = i + 4Z + 4Z* + 2Z 3 , C 2 (Z) = 2 + 2Z + 3Z 2 + 2Z 3 , and C 3 (Z) = 0. Hence c 0 0 = 1, c 01 = 2, c Q 2 = 1, c 0>3 = O,c 10 = 2, 

C 1 ,1 ~ ^’ C l,2 ~ ^’ C l,3 = 2 > C 2,0 = 2,C 2,1 = 2,C 2,2 = 3 ’ C 2,3 = 2 > C 3,0 “ C 3,l = C 3,2 = 3I1C * C 3,3 = 


100 



o 
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